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We present a microscopic theory of transport through quantum dot set-ups coupled to super- 
conducting leads. We derive a master equation for the reduced density matrix to lowest order in 
the tunneling Hamiltonian and focus on quasiparticle tunneling. For high enough temperatures 
transport occurs in the subgap region due to thermally excited quasiparticles, which can be used 
to observe excited states of the system for low bias voltages. On the example of a double quantum 
dot we show how subgap transport spectroscopy can be done. Moreover, we use the single level 
quantum dot coupled to a normal and a superconducting lead to give a possible explanation for the 
subgap features observed in the experiments of Ref . 1 . 
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I. INTRODUCTION 

In the last two decades modern fabrication techniques 
made it possible to connect quantum dot systems with 
superconducting leads. Quantum dots were realized 
with carbon nanotubes 1 7 , metallic particles 8 , semicon- 
ducting nanowires 9 12 , single fullerene molecules 13 , self- 
assembled nanocrystals 14 and graphene quantum dots 15 . 
The experiments show a gap in the Coulomb diamonds 
which is proportional to the superconducting gap, reflect- 
ing the BCS-density of states. In the sequential tunneling 
regime higher order quasiparticle tunneling processes are 
suppressed and current flows due to single quasiparticle 
tunneling. First transport theories were presented 16 , us- 
ing a master equation approach, where the rates were 
calculated on the basis of Fermi's golden-rule. Another 
method based on non-equilibrium Green's function was 
used by Yeyati et al. 17 and K. Kang 18 to describe res- 
onant tunneling through an effective single level quan- 
tum dot in the limit of very strong Coulomb repulsion 
in the dot (U — >• oo limit), where transport is governed 
by quasiparticle tunneling; the corresponding I-V curves 
show an intrinsic broadening of the BCS-like feature in 
the current in agreement with experimental observation 8 . 
For small Coulomb repulsion, higher order processes lead 
to Josephson current 9 and Andreev reflections 2 5,7,10,15 , 
which appear as subgap features in the experiments. 
Both effects were studied intensely experimentally and 
theoretically 4 ' 17, 19 ' 20 and were recently summarized in re- 
view articles of Refs. 21 and 22. Besides Andreev reflec- 
tions also the Kondo effect 13 as well as Yu-Shiba-Rusinov 
bound states 5,23,24 can lead to subgap features and are 
subject of current research. If the temperature becomes 
comparable with the superconducting gap quasiparticles 
can get thermally excited across the gap, leading to ad- 
ditional subgap features 16 . 

In the following we present a microscopic theory for 
transport through superconducting hybrid nano junctions 
for finite superconducting gap |A| < oo in the sequential 
tunneling limit. In particular, we trace out all degrees 
of freedom of the superconducting leads to obtain a gen- 



eralized master equation for the reduced density matrix 
to lowest order in the tunneling Hamiltonian. In con- 
trast to Green's function techniques, see e.g. Ref. 22, 
this method enables to treat the interactions on the sys- 
tem exactly. Moreover, as shown on the example of a 
double quantum dot, our theory can treat any quan- 
tum dot set-up if its many-body energy spectrum and 
eigenstates are known. We focus on transport involving 
thermally excited quasiparticles, and show that excited 
states of the quantum dot system can be observed in the 
current voltage spectroscopy in the Coulomb blockade re- 
gion. Though transitions between two ground states are 
blocked due to the gap in the BCS-density of states, ther- 
mally excited quasiparticles can participate in transport 
through excited system states, giving a source of subgap 
features in superconducting hybrid systems. For a quan- 
tum dot coupled to a normal and a superconducting lead, 
a possible explanation for the subgap features observed 
in Ref. 1 is given, where a carbon nanotube quantum dot 
is coupled to a normal and a superconducting contact. 

The manuscript is organized as follows: In Sect. II we 
introduce the Hamiltonian in a system-bath model using 
a number conserving version of the Bogoliubov-Valatin 
transformation 25,26 to describe the electrons of the su- 
perconducting leads as quasiparticle excitations of the 
BCS-ground state and Cooper pairs, introducing Cooper 
pair creation and annihilation operators. The Cooper 
pair operators are essential to conserve the total particle 
number and have important consequences on the trans- 
port properties. In Sect. Ill we derive an explicit form 
of the Cooper pair operators following Ref. 27 and dis- 
cuss the influence on thermodynamic properties of the 
leads. In Sect. IV, the general master equation for the 
reduced density matrix is derived and used to calculate 
the current. In Sect. V we use the theory to calculate 
transport characteristics for two systems, the single level 
quantum dot (SD) and the double quantum dot (DD), 
the later in two possible configurations cf. Fig. 1. The 
SD is used to explain basic phenomena like a gap open- 
ing in the Coulomb diamonds which is proportional to 
the superconducting gap, and transport involving ther- 
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Figure 1. Sketch of the transport set-up of a double quantum 
dot (DD) coupled to superconducting leads. The DD is illus- 
trated in the parallel (top panel) and serial (bottom panel) 
configuration. Tunneling events are depicted by arrows. 



mally excited quasiparticles 16 . On the other hand, the 
DD possess a richer many-body spectrum with several 
excited states. We visualize transitions through excited 
system states in the low bias regime using thermally ex- 
cited quasiparticles. Due to the gap in the BCS-density of 
states, the ground state to ground state transition is not 
allowed in all cases, leading to transport through excited 
system states, appearing as peaks in the Coulomb block- 
ade region. The threshold for observing excited system 
states in the subgap region is that the energy difference 
between the excited state and its ground state must be 
smaller than 2|A|. We confirmed this threshold by means 
of the independently gated DD, where the detuning of the 
two sites changes the level spacing. Finally the N-QD-S 
system is investigated, where a quantum dot is coupled 
to a normal and a superconducting lead. In this case only 
the superconducting lead produces thermal lines in the 
Coulomb blockade region, giving a possible explanation 
for the subgap features in Ref. 1. 



II. MODEL HAMILTONIAN 

In the following we consider quantum dot systems 
weakly coupled to two superconducting leads. The to- 
tal Hamiltonian is written in a system-bath model: 



H = H. 



(i) 



where Hs represents the Hamiltonian of the quantum dot 
system, Hb is the Hamiltonian of the superconducting 
leads, and Ht describes the tunneling between the sys- 
tem and the leads. Specifically, we focus on two systems, 
a single level quantum dot (SD) and a double quantum 
dot (DD). The SD has been focus of many theoretical 
works before 16 20 , and we use its simple Fock-space struc- 
ture to demonstrate some generic effects resulting from 
the superconducting leads. 



We describe the SD by the single impurity Anderson 
model: 



H S d = ^2e d dl d a +17 



(2) 



where n a = d a d a is the number operator of the elec- 
trons on the dot with spin cr. This model describes a 
quantum dot with on-site energy and Coulomb repul- 
sion U which can be occupied by at most two electrons. 

The highest occupied state is defined as |2) = d^ d^ |0), 

the one particle states are defined as \la) — d^ |0), and 
|0) is the state with zero particles. 

For the DD we use a modified version of the Pariser- 
Parr-Pople Hamiltonian 28,29 : 



H D D = ^2 6a(T ^a* + ^2 ( b ^ <W +&* d la 

a £{1,2} cr ^ 



(3) 



Here, d ao . are the creation operators for an electron on 
site a G {1,2} with spin cr. They define the number op- 
erators h acr = d} aa d aa . The operator n a = h a ^ + n a ± 
counts the number of electrons on site a. In the gen- 
eral case we distinguish between the four on-site energies 
e aa and between the on-site Coulomb interactions U a . 
Electrons on different sites interact through the inter-dot 
Coulomb interaction V; b describes the hopping between 
the two sites. In our set-up the on-site energies can be 
controlled by capacitively coupled gate electrodes. In 
the case of site-independent on-site energies and on-site 
Coulomb interaction the Hamiltonian can be diagonal- 
ized analytically 30 ' 31 . 

The superconducting leads are described by the mean 
field form, H B dF of the pairing Hamiltonian, where we 
additionally inserted a unity represented by a product of 

Cooper pair annihilation and creation operators, S^ S^ = 
1, which will be specified later in Sec. III. We find 



rjka 



^ (A?? cj^ C^_ k ^ S^ +A^ S^ Crj- 



rjk 



(4) 



where ^ = e& — measures single particle en- 
ergies ek with respect to the electrochemical poten- 
tial /i^, and N v = ^2 kcr c^ k(T c v k a counts the number 
of electrons in lead rj. Finally, A^ = lA^e 2 ^ = 
- IZi Vik $1 C77-/4 c,jfct) denotes the superconducting gap 
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of lead rj. Here (•) denotes a thermal average calculated 
self-consistently using the mean field Hamiltonian of Eq. 
4. 

The tunneling Hamiltonian, 

Ht = ^2 tr ^ aa ^lk* da* "f^a daa <W> (5) 
rjk aa, 

describes the tunneling between the leads and the two 
sites of the DD, where the tunneling coefficients t vcca de- 
pend on the lead, site, and spin index. Depending on the 
choice of the tunneling coefficients the DD is described 
in parallel or in serial configuration, see Fig. 1. For the 
single dot we skip the index a in Eq. (5), as only one site 
is involved. 



A. Diagonalization of the lead Hamiltonian 



where <p v is the phase of the superconducting gap A^. 
Applying the transformation of Eq. (6) on Eq. (4) we 
obtain: 

H B = Y J E vk l\ ka l vka +E G + Y,^K (12) 

rjka 7] 



In Eq. (12) E^ = y + 2 denotes the quasipar- 

ticle energy, and Eq is a constant energy off- set, often 
referred to as the energy of the Cooper pair condensate. 
For later reference we note that the term ^2 is 
not included in the diagonalization procedure and is still 
written in terms of electron operators. 

III. THE BCS GROUND STATE AND ITS 
CONNECTION WITH THE COOPER PAIR 
OPERATORS 



The most famous way to diagonalize the mean field 
Hamiltonian, Hq, of Eq. (4) was first introduced 
by Bogoliubov 32 . We are following Josephson and 
Bardeen 25,26 who modified the so called Bogoliubov 
transformation in a number conserving way. We adopt 
this idea and define the Bogoliubov transformation: 



(6) 



where a = —a. In Eq. (6) creates a fermionic 

quasiparticle, often called bogoliubon, which is defined 

by 



% ka |GS>„ = 0. 



(7) 



(8) 



Here |GS) denotes the BCS ground state, or Cooper pair 
condensate of lead rj. Bogoliubons are quasiparticle exci- 
tations of the Cooper pair condensate, meaning that the 
Cooper pair condensate is defined as the vacuum state of 
the bogoliubons, see Eq. (8). The Cooper pair creation 

operator S^ creates a Cooper pair in lead 77. It has to be 
defined in such a way that the unitarity of the Bogoli- 
ubov transformation is conserved. The coefficients u v k 
and are complex numbers and fulfill: 



1. 



They read: 




v v k 




Cqk 

\E v k\ 



(9) 



(10) 



(11) 



In the microscopic description of superconductive tun- 
neling it is necessary to know the analytical form of the 
Cooper pair operators. In this section we show the con- 
nection between the Cooper pair operators and the BCS 
ground state by rewriting the BCS ground state as a 
coherent superposition of states with a fixed number of 
Cooper pairs. We start from 33 : 



|gs) = JJK + ^r1)|o), 



f - At At 



C kt C -ki' 



(13) 



Ex- 



where |0) is the vacuum state and R 
ploiting the fermionic properties of the electron operators 
we can write: 

|GS) = n^exp(^^Ri) |0> 

h \ u U k J 

(14) 

Eq. (14) can be written in a more compact form 27 , 

00 

\GS) = J2 b n\n), (15) 

by defining 



n=0 



10} • 



(16) 



(17) 



In App. A 1 we prove that \n) contains 2n electrons which 
form n Cooper pairs. In Eq. (17) we have introduced an 
additional coefficient 



E 



M 
Wk\ 



n/2 



(18) 
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to normalize the states: 



where the trace is over many-body quasiparticle states 



(n\m) = S n 



(19) 



Following Eq. (17) we define the Cooper pair creation 
operator: 



(20) 



such that a state with n Cooper pairs is created by ap- 
plying the Cooper pair creation operator n times on the 
vacuum: 



|n) = (S T ) n |0). 



(21) 



As we see from Eqs. (15) and (21) the ground state of the 
bogoliubons is created by Cooper pair creation operators. 
Using Eq. (20) we can proof the following properties of 
the Cooper pair operators: 



and 



[S',S] = 0, 



S S = 1, 



(22) 



(23) 



(24) 



(25) 



From Eq. (25) it follows that the bogoliubon operators 
and the Cooper pair operators commute. 



A. Thermodynamic properties of the leads 

The description of electrons in terms of bogoliubons 
and Cooper pairs makes it necessary to discuss the ther- 
modynamic properties of the superconducting leads. The 
chemical potential in this description is the same as for 
the electrons, since the unitary transformation conserves 
all electronic degrees of freedom. In the following we drop 
the lead index 77, and consider only one lead. 

In order to calculate thermal expectation values we use 
the equilibrium density matrix of a superconductor: 



PB 



Z G 



(26) 



where Hq = Hb — pN , f3 is the inverse thermal energy, 
and Zq is the partition function in the grand canonical 
ensemble. We find that the thermal expectation value of 
a pair of Bogoliubov quasiparticles is equal to the Fermi 
function: 



Tr B (tL Ik* PB 



1 



e (3E k + I 



~m, (27) 



\{n}) = \n kiai ,n k2a2 ,...)=A ]J tf qr |GS) . (28) 

(qr)e{n} 

The factor A normalizes the many body state of Eq. (28) , 
and {n} is a set of occupation numbers. We can see that 
tracing over the many body states includes the Cooper 
pair degrees of freedom through the ground state. Hence 
the trace includes all electronic degrees of freedom. The 
derivation can be found in App. B. 

For later reference, it is important to know thermal 
expectation values of combinations of bogoliubon and 
Cooper pair operators, for instance: 



TrB(S f 7L7 fe0 



p B )=f + (E k ). 



(29) 



In fact, since the Cooper pair operator is commuting with 
the bogoliubon operators, see Eq. (25), we need to calcu- 
late the BCS-ground state expectation value of a single 
Cooper pair operator, which is equal to one in the ther- 
modynamic limit, see App. B: 



(GS|S |GS) = 1. 



(30) 



IV. TRANSPORT THEORY AND THE 
GENERALIZED MASTER EQUATION 

The expectation value O = (O) = Tr (Op) of any ob- 
servable associated to an operator O can be evaluated 
once the total density operator p is known, cf. Eq. (46) 
below. To this extent we start from the Liouville equa- 
tion for the density operator in the interaction picture, 
see e.g 



34. 



,9 . , N 



(31) 



Eq. (31) can be formally integrated and reinserted back 
into itself, 



ih'p^t) = [H T j(t)Ji(0)] 



H T At),[HT,i(t'),Pi(t') 



(32) 



which is still exact at this level and allows a perturbative 
treatment in the tunneling Hamiltonian Ht- 

Prior to time t = the bath and the system do not 
interact, meaning that the total density matrix of system 
and leads are decoupled: 



Pi(0)=Ps(0)pb(0). 



(33) 



The density matrix of the leads, can be described by 
the equilibrium thermodynamic expression shown in Eq. 
(26). Further we assume that the leads have so many 
degrees of freedom that they stay in thermal equilibrium 
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up to a correction of order H T . It is convenient to trace 
out the degrees of freedom of the leads and define the 
reduced density matrix: 



Predjif) EE Tr B pl(t). 



(34) 



In the Schrodinger picture, the master equation for the 
reduced density matrix reads: 



Pred{t) = -r[pred(t),H S ] 



U (t) / dt'x 



x Tr B 



H T ,i(t), 



H T j(t),pred,l(t')p L 



(35) 



where we neglect terms of order O(H^) and XJo(t) = 
e -jrH s t - g ^Yiq time evolution operator of the unperturbed 
system. 



A. Superconducting leads 

The features of the superconducting leads are revealed 
when using the Bogoliubov transformation (6) to express 
the tunneling Hamiltonian. This yields additional terms 
compared to normal conducting theory. 



1. Time evolution of the quasiparticles 

To proceed we have to specify the time evolution of 
the Bogoliubov and Cooper pair operators. We find: 



Inka' 



(36) 



(37) 



in agreement with the results of Josephson and 
Bardeen 25 ' 26 . When calculating the time evolution it is 
important to remember that in the lead Hamiltonian the 
term p^N^ is still written in terms of electron operators. 

Before we proceed, we like to emphasize the impor- 
tance of the Cooper pair contribution for finite bias volt- 
ages. As already pointed out by Governale et a/. 20 , in 
this case Pt\ 

cannot be set to zero and the time evolution 
of the Cooper pair operators, Eq. (37), plays an impor- 
tant role. Neglecting the Cooper pair contribution for 
finite bias voltages 35 violates the number conservation in 
the tunneling processes and can lead to coherences which 
would vanish in the number conserving case. 



2. Difference to the normal conducting theory 

To compute Eq. (35) we rewrite the electron operators 
using the Bogoliubov transformation, Eq. (6), and insert 



the time evolution as in Eqs. (36) and (37). This yields 
four different traces to be calculated. We find: 

( 38 ) 



Tr B l^nkajit) cl, k , a , j(t')p B 

(39) 

+ |^ fc |V + (^fc)e + * ( ^ fc -^ )(t - t,) |, 

sgna^^e^ 2 ^ £ |/-(^)e-^^+^)^-^) (40) 
- /+(£; r7/c )e + ^ ( ^ fc "^ )(t "^ ) 



sgna ^^ fc e-^ 2 ^ t |/ + (^)e+^^ + ^)^- t/ ) (41) 



/ (E v k)e 



-jr(E vk -H v )(t-t 



"}■ 



where f~(E) = 1 — f + (E). In contrast to normal con- 
ducting theory, the traces of Eqs. (40) and (41) are not 
vanishing in the superconducting case. However, to sec- 
ond order in perturbation theory they do not contribute 
to the time evolution of the density matrix, for reasons 
to be presented later. They are crucial in higher order 
perturbation theory and yield e.g. the phase-dependent 
Josephson current 20 . 



B. General Master Equation for the reduced 
density matrix 

Collecting all the previous results and expressing Eq. 
(35) in the basis of the system eigenstates we obtain 
the Bloch-Redfield form of the general master equation 
(GME) for the reduced density matrix: 

Pnn' = — ^ {E n — E n >)p nn '(t) 



El , r>N^N-l 

\ nn'mm' ' nn'mm' 



(42) 



nn'mm' * ^ nn'mm' \Pmm'\ L )i 
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where we defined the Redfield- tensors: 



nn' mm' / ^ 



5m' ri ^ {^nllm)r] ^ mn fim'lln' 



N^N±1 
V 



\ m'n'nmJrj \ m'n'nmJrj [' 



(43) 



and 



^nn'mm' / v 



iV— >-iV±l 
'77 



V m'n'nm)r] V m'n'nm/rj (' 



(44) 



In Eqs. (43) and (44) we introduced rates T, which 
we call "normal", and rates E, which we call "anoma- 
lous" . The anomalous rates E come from terms contain- 
ing traces of the type of Eqs. (40) and (41). In App. C 
we show that they are equal to zero. The normal rates 
originate from terms containing traces of the type of Eqs. 
(38) and (39). Further, we distinguish between rates de- 
scribing the increase and rates describing the decrease of 
the particle number on the system, emphasized with the 
superscript N —> N± 1. Their detailed form is presented 
in App. C. The normal rates with the superscripts ± are 
connected by complex conjugation and reversing of the 
indices: 



,N^N±l 



(r + , 

V n'r 



(45) 



C. Current 

Having derived the GME for the reduced density ma- 
trix in Eq. (42), we can use it to calculate measurable 
quantities like the current and the differential conduc- 
tance. In this section we present an expression for the 
current derived from the second order GME of Eq. (42). 
To do this we introduce a current operator whose statis- 
tical average gives the total current: 



Irj = Tr (Irjpto 



(46) 



In general, the current operator of lead 77 is defined as 
the variation of the total particle number in lead 77 with 
time: 



-ejN v At) 



+ie 



N vJ (t),H TJ (t) 



(47) 



Calculating the commutator of Eq. (47), we see that the 
current operator has the same operatorial structure as 
the tunneling Hamiltonian: 



ha. 



(48) 



differing only in the prefactor and summation. Hence, 
by applying the same perturbation theory as before, we 
obtain for the current in lead 77: 



V*) = eE((rr +1 ) r (r 

nml 



N^N-1\ \ N 



/£»(*)■ (49) 



In Eq. (49) we introduced the abbreviations 

/ piV— >-iV±l\ = /p+ ^N^N±l / ,N^N±1 

\ nmm'n' ) rj \ nmm'n' J rj V nmm'n' J rj 



= 2Re((r+ Wn ; 



N—*N-\-l 

TJ 



(50) 



exploiting Eq. (45). This gives us rates which are real 
and read: 



(C^ 1 ), = Re[ K mm , n , D(E m , n , - N + i 7 ) 



(51) 



(r 



nmm'n' ) 



v = Re\Tl, n , nm D{E n>m ,-n v + i^ 
x f~ [E n i 



(52) 



where 



2tt 



Kmm'n' = W^a'a ( n l \™) (m'\ d a , a \ri) . 

aaa' 

(53) 

In Eqs. (51) and (52) E n r m > = E' n —E' m denote differences 
between system eigenenergies and 



D(E) = p N Re 



\E\ 



(54) 



is the BCS-density of states, with pN = 27^2 labeling 
the density of states for normal leads which is assumed 
to be constant around the Fermi level; V denotes the vol- 
ume of the lead and m is the electron mass. In order to 
renormalize the divergence of the density of states we in- 
troduced a finite lifetime h/j of the quasiparticle states in 
the superconducting leads, leading to a Lorentzian broad- 
ening of the resonance condition, see App. C 3. This 
assumption is also in agreement with the results of Levy 
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Yeyati et al 17 , where they showed that the broadening of 
the BCS-like features in the current is due to the coupling 
to the leads. Eq. (49) is a general result and can be ap- 
plied to any transport set-up where an arbitrary system 
with discrete levels is weakly coupled to superconducting 
or normal conducting leads. The normal conducting case 
is obtained by setting |A^| = and 7 = 0. 

In this article we are only interested in the stationary 
limit. Hence, we replace the density matrix in Eq. (49) 
by its stationary solution which is determined from Eq. 
(42) by imposing p^ n , = 0. 



V. TRANSPORT THROUGH MULTIPLE 
QUANTUM DOT DEVICES 

In the preceding sections we developed a perturbative 
microscopic theory for the stationary current of quantum 
dot devices coupled to superconducting leads. In the 
following, we show the predictions of the theory for two 
models, the single level quantum dot (SD) and the double 
quantum dot (DD). In the transport set-up the bias and 
gate voltages influence the energy configuration of the 
leads and the system, respectively. Specifically, the bias 
voltage is modifying the electrochemical potential of the 
leads, which we choose to have a symmetric voltage drop. 
Therefore we define the chemical potentials of the left and 
right lead, respectively: 



Vl/r = Mo ± e 



V b 



(55) 



where /xq is the equilibrium chemical potential. The gate 
voltages are modifying the on site energies of the system: 
we replace — > + eV g in the SD- and e a —> e a + eV g a 
in the DD-Hamiltonian. Here e = — \e\ is the electron 
charge. 

In the following we neglect coherences in the GME, 
considering only diagonal contributions of the reduced 
density matrix p nn by setting n = n' in Eq. (42). Hence, 
it suffices to use only two indices for the transition rates. 

In current voltage spectroscopy it is convenient to il- 
lustrate the conditions under which current is allowed 
to flow as lines in the stability diagrams. These so called 
transition lines are fixed by the energetic part of the tran- 
sition rates at the source 77 = S and the drain rj = D 
contact: 



(r 



N—^N-\-l 



) <xf+(AE-^)D(AE-n v ), (56) 
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Figure 2. Panels (a) and (b): Density of states (continuous 
line) and Fermi function (dotted line) at ksT = 0.2 meV and 
ksT = 0.01 meV, respectively. Panels (c) and (d): Product 
of the density of states and the Fermi function for the tem- 
peratures used in Fig. (a) and (b), respectively. 



eV b 



AE n = AE„ > \A\ 



Dt-D 




-iV+J 



st-s 



D+Dt+ 



Figure 3. (Color online) Illustration of the transition lines 
appearing in presence of superconducting leads. The green 
lines mark transitions at the Source and the Drain contacts, 
described by the inequalities of Eqs. (58), (59), (62), and (63). 
The red lines mark transitions involving thermally excited 
quasiparticles, given by Eqs. (60), (61), (64), and (65). The 
Eg-N diagrams for the points (a)-(e) are sketched in Fig. 5. 



(r^ 1 ^),, oc f~(AE - [i v )D(AE - (i v ), (57) 

neglecting the lifetime broadening 7 for simplicity, and 
with AE = E^ +1 — E^ the energy difference of the two 
transport levels. Fig. 2 illustrates this product for two 
different temperatures: for high enough temperatures 
quasiparticles can be excited thermally across the gap 
giving a small peak in the transition rates 16 . The peak 



positions define transition lines when plotted in a V g -Vb 
diagram. Notice that while the most pronounced peak 
survives also at zero temperature and defines a transport 
threshold, the second peak vanishes at low temperatures 
and essentially only processes at and close to the peak 
are relevant. For an N — >> N + 1 transition we denote 
transitions associated to the more pronounced peak as 
S+ and D+ when happening at the source or at the drain 
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contact, respectively. Transitions involving thermally ex- 
cited quasiparticles are called St+ and Dt+. In complete 
analogy, we classify transitions from N + 1 — >• N: we de- 
note by S- and D- the more pronounced transitions at the 
source and at the drain, and by St- and Dt- their thermal 
counterparts. In total we find 8 different transition lines, 
as depicted in Fig. 3. In the following we derive trans- 
port conditions and provide equations for the transport 
lines. For convenience we introduce AE g = AE — /i . 

We start with the analysis of the N —> N + 1 transi- 
tions, which are described by the rates in Eq. (56). From 
the arguments we find that the rates do not vanish if 

Source S+: 



AE g < -\A\ 



Drain D+: 



AE a < 



eVb 
2 



eVb 
2 



(58) 



eV b 



< -AE g - \A\. (59) 



Another contribution comes from the thermally excited 
quasiparticles states, namely, if the argument of the 
Fermi function f + (AE — n v ) and of the density of states 
D(AE — fi v ) is equal to |A|. At this point the transition 
rates are peaked and contribute to the current: 
Source (thermal) St+: 



AE g = \A\ 



en 

2 



eVb 
2 



= AE n 



Drain (thermal) Dt+: 



AE a 



eVb 
2 



2 



-AE„ + |A|. 



(60) 



(61) 



Since the thermally excited quasiparticles produce a peak 
rather than a step in the current voltage characteristic, 
the corresponding transport condition is formulated with 
an equality. 

Transitions from N + 1 — >> TV are described by the rate 
of Eq. (57), leading in complete analogy to the previous 
case to the following transport conditions: 

Source S-: 



-AE g < 
Drain D- 



-AE g < 



eVb 
2 



2 



eV b 



< AE n 



eV h 



> -AE„ 



Source (thermal) St- 



-AE„ 



eV b 



■ e -^=AE g + \A\ 



(62) 



(63) 



(64) 



Drain (thermal) Dt- 



eV b = 



eV b > |A| 



| Source 
| Drain 



. St+ 
; Dt+ 



% Source 
1 Drain 



eV b < |A| 



eV b = |A| 



. st+ 

- Dt+ p 



| Source 
| Drain 



. St+ 
. Dt+ 



^ Source 
I Drain 



N-l N N+l 



eV b = 2IAI 



N-l N N+l 



eV b > 2IAI 



■. Dt+ m s ° urce 

^ Drain 



^ Source 
1 DL+ ^ Drain 



Figure 4. (Color online) Visualization of the transport condi- 
tions of Eqs. (58)- (65). We plotted the threshold of the trans- 
port inequalities as green lines (S=b, D±); for the equalities 
coming from transitions involving thermally excited quasipar- 
ticles we used red lines (St=t, Dt± ). Choosing the reference 
level in the N particle subspace, we found a scheme where 
transitions are energetically allowed to levels which lie in the 
shaded region below the green lines and to levels which align 
with the red lines. The dashed box marks the bias window 
eV b . 



1 . Visualization of the transport conditions 

To visualize the transport conditions of Eqs. (58)- (65) 
we extend the scheme of Donarini et al. of Ref. 36 to 
superconducting leads. The scheme is depicted in Fig. 
4 and illustrates for which relative position of the sys- 
tems eigenenergies = E^ — /iqN transitions are en- 
ergetically allowed. The green lines mark the borders of 
the inequalities, the red lines the sharp equalities for the 
thermal transitions, meaning that transitions can occur 
to states lying below the green lines, and to states which 
coincide with the red lines. In order to see a transition 
between two levels in the stability diagram a source and 
a drain transition must be allowed between the two levels 
(depicted as arrows in the E g -N diagrams of Fig. 5). We 
note that for a full analysis of the transport properties 
also the geometrical part of the rates must be taken into 
account and transport occurs only if T ^ 0. 



A. Single level quantum dot model 



-AE g = \A\ 



eVb 
2 



eVb 
2 



-AE g - \A\. 



(65) 



The simplest example of a quantum dot system is the 
single level quantum dot presented in Eq. (2). Since only 
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Figure 5. (Color online) (a)-(d): E g -N diagrams for a single 
level quantum dot with AE g > |A| and at bias voltages as 
sketched in Fig. 3. For the simulations of Fig. 6 AE g > A 
corresponds to a gate voltage eV g < — 2.6meV. In (a) we cut 
the S+ line: the particle number on the system is increased by 
a tunneling event at the source contact and decreased at the 
drain, (b) Cut with the thermal line St+: the particle num- 
ber of the system is increased by a tunneling event involving 
a thermally excited quasiparticle at the source contact and 
decreased by tunneling into empty states in the source drain, 
respectively. Increasing further the bias voltage until its sign 
changes leads to an interchange of the role of the source and 
the drain voltage. Hence, (c) and (d) describe the same sit- 
uation as (a) and (b) under S D. (e): E g -N diagrams for 
a single level with < AE g < |A|. The two levels are only 
connected by two drain transitions, meaning that in this con- 
figuration the system is in thermal equilibrium with the drain 
contact. 



one level is involved, we can do most calculations analyti- 
cally and understand the basic mechanism resulting from 
the superconducting leads. In Fig. 6 the stationary cur- 
rent is shown as a function of bias and gate voltage for 
superconducting leads at k B T = 0.5|A|. We observe the 
expected gap 5 between the Coulomb diamonds which is 
equal to 4|A|/e. The gap can be explained using Fig. 
3 and the corresponding Eqs. (58)- (65). One dashed 
line marks the gate voltage where AE g = 0. Along 
this line the conditions under which current is allowed 
to flow read: eVb/2 > |A| for the S+, D- lines, and 
eVfc/2 < — | A | for the S-, D+ lines, opening a bias window 
of 4|A|/e where current is blocked for low temperatures 
k B T <C |A|. For higher temperatures of k B T w 0.5|A| 
we observe small peaks in the Coulomb blockade region 
(green area) which are due to thermally excited quasi- 
particles, they correspond to the red lines in Fig. 3. In 
Fig. 5 we show the energy particle number diagrams in 
the points (a)-(d), which lie on a vertical cut through 
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Figure 6. (Color online) (a) Current voltage characteristics 
of a SD coupled to superconducting leads. Parameters are 
ksT — 0.3 eV and |A| = 0.6 meV. (b) Subgap features com- 
ing from thermally excited quasiparticles of the 0-1-particle 
transition, highlighted as a dashed box in (a). 



Fig. 3 at AEg > \A\ which corresponds to a gate volt- 
age eV g > 2.6 meV in Fig. 6. In Fig 5 (a) we depicted 
the Eg-N diagram for a cut with the S+ resonance line, 
where the particle number on the system is increased 
by a tunneling event at the source and decreased at the 
drain contact. For bias voltages smaller than the one at 
resonance (corresponding to larger eV& as e is the nega- 
tive charge of an electron) the S+, D- transitions remain 
open and current can flow. In Fig. (5) (b) the E g -N dia- 
gram at the resonance line St+ is shown. In this case the 
bias voltage is not large enough to allow the transitions 
S+ of Eq. (58). For low temperatures no quasi par- 
ticle is thermally excited meaning that only transitions 
from 1 — >• are energetically allowed (green arrows). For 
high enough temperatures, however, the particle number 
of the system can be increased by tunneling events in- 
volving thermally excited quasiparticles opening the St+ 
transition. By changing the sign of the bias voltage the 
role of the source and the drain is inverted, explaining 
the transition lines Dt+ and Dt (Fig. 5 (c) and (d)). 

Another interesting constellation of the energy level 
occurs in the region of < AE g < |A| (Fig. 5 (e)), 
where in the current-voltage characteristics the thermal 
lines are vanishing. Transitions can only occur at the 
drain contact, as the bias is not large enough to allow 
transitions at the source. Hence, the system is in thermal 
equilibrium with the drain contact and the occupation 
probabilities are related by the Boltzmann distribution: 



Po _ 
Pi 

in the limit of 7 — >• 0. 



o 0(AE g +eV b /2) 



(66) 



B. The double quantum dot 

We have seen that the theory can reproduce well known 
results for the SD and we understood the properties of 
the thermal transitions in E g -N diagrams with only one 
non degenerated level per particle number. In the fol- 
lowing we investigate a more advanced system, the dou- 
ble quantum dot, where the many body spectrum gives 
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rise to more than one non degenerated level per parti- 
cle number, so called excited system states. For normal 
conducting leads the excitations cannot be seen for low 
bias voltages, since transitions to the ground state are 
always possible, blocking transport through the excita- 
tions. In the last subsection we have seen that for super- 
conducting leads the energy difference must be at least 
\AE g I > eVb/2— | A | to have non thermal source and drain 
transitions. Hence, we find situations where the transi- 
tion to the ground state are energetically not allowed and 
transport occurs through excited system states. 

We start with equally gated dots with the same on- 
site energies and on-site Coulomb interactions, where it is 
possible to diagonalize the Hamiltonian analytically 30,31 . 
In the second part, the case of independently coupled 
dots is discussed, where the detuning of the two gate volt- 
ages influences the level spacing of the energy spectrum. 
Thus, excited states can be observed only in detuning 
ranges where the difference between the energy level of 
the excited state and its ground state is less than 2|A|. 



1. Equally gated dots 

For equally gated dots the on-site energies of the two 
sites are modulated with the same gate voltage. Hence, 
it is convenient to plot the current as a function of the 
bias and the gate voltage as for the SD. Fig. 7 shows the 
current of an equally gated DD in serial configuration. 
As for the SD we observe Coulomb blockade and the gap 
of 4|A|/e between the tips of the diamonds. Transport 
carried by thermally excited quasiparticles is of partic- 
ular interest, as it allows to observe transitions through 
excited system states for low bias voltages, which are 
often diminished by the ground state transitions in the 
normal conducting case. In order to show some interest- 
ing phenomena resulting from the more complex spec- 
trum, we concentrate on the 0- to 1-particle transition 
where three levels are involved. In the 1-particle spec- 
trum, the difference between the ground state and the 
excited state is equal to 2|6|, where b < is the tun- 
neling strength between the two dots. Meaning that by 
tuning the coupling between the two dots it is possible 
to influence the level spacing. Fig. 8 shows a sketch of 
the transition lines expected for the — 1 transition for 
\b\ < |A|, where the red (green) lines show the ground 
state to ground state transitions, and the blue (orange) 
lines the ground state to first excited state transitions. 
For a better understanding of the transport properties we 
cut the transitions lines horizontally for a small bias volt- 
age eVb/2 < |A| in the Coulomb blockade region (points 
(A)-(D)), the corresponding E g -N diagrams are depicted 
in Fig. 9. In point (A) the difference between the ground 
states is equal to AE g = eVb/2 + |A| opening the ther- 
mal transition St+ and current can flow. Following the 
dashed line to point (B), the 1 particle states are shifted 
down in energy until the St+ transition is allowed be- 
tween the 0-particle ground state and the 1-particle ex- 
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Figure 7. (Color online) (a) Current voltage characteristics of 
an equally gated DD in serial configuration at ksT = 0.2 meV 
and |A| = 0.4 meV. (b) I-V characteristics in the subgap 
region corresponding to the dashed box in (a). The distance 
between the 1-particle excited state and its ground state is 
equal to the coupling strength 2\b\ of the two dots. Moreover, 
2\b\ < 2|A|. The black arrow marks the transition line coming 
from transport through the 1-particle excited state, (c) I- 
V-characteristics in the subgap region, where we increased 
the coupling between the two dots, leading to a level spacing 
which is larger than 2|A|, hence transport through the excited 
system state is not allowed and the line disappears. 



cited state. Since |6| < A, the 1-particle ground state is 
energetically not accessible and current can flow through 
the excited state. We like to emphasize that the blocking 
of the ground state transition is only valid as long as the 
distance between the two 1-particle levels is smaller than 
2|A|. For larger distances the ground state is energeti- 
cally accessible, blocking the current through the excited 
state, c.f. Fig. 10. In point (C) eV g is further decreased, 
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Figure 8. (Color online) Sketch of the transition lines for the 
0-1 particle transition of an equally gated DD. It shows two 
copies of Fig. 3 where the labeling of the blue (orange lines) 
is the same as for the green (red) lines. The blue (orange) 
lines mark the transition lines corresponding to the 0- particle 
ground state to 1-particle first excited state transition. 
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Figure 9. (Color online) E g -N diagram corresponding to the 
points of Fig. 8 where the dashed line cuts the transition lines 
for the case of an equally gated DD. In this case the distance 
between the 1-particle ground state to the 1-particle first ex- 
cited state is equal to 2b < 2|A|, where b is the tunneling 
strength between the two quantum dots. (A) Point on the 
thermal line St+ of the ground state to ground state transi- 
tion. (B) Point on the thermal line St+ of the ground state to 
first excited state transition. (C) Point on the Dt- line of the 
ground state to ground state transition. (D) Point on the Dt- 
line of the ground state to first excited state transition; this 
line cannot be seen in the current voltage characteristics, as 
the ground state to ground state transitions are open. Hence, 
in the long time behavior the system will occupy the 1-particle 
ground state blocking the current through the excited state. 



st+ CB) 

_[D \|jeV»/2 + |A| 
> 



1 N 

Figure 10. (Color online) E g -N diagram of point (B) in Fig. 
8, for a level spacing of the one particle energies larger than 
2b > 2 A. In contrast to Fig. 9 the transition between the 0- 
particle ground state and the 1-particle excited state is open, 
blocking the current. 



the Dt- transition between the ground states is opening, 
and current can flow. Point (D) shows the typical energy 
configuration in which current through the excited state 
is blocked, even though the transition through the ex- 
cited state is energetically allowed. The reason for that 
is the 1-particle ground state which can be populated, 
but transitions describing its depopulation are energeti- 
cally not allowed, leading to a blocking of the current in 
the stationary limit. 

To demonstrate the important role of the level spacing 
we show the current voltage characteristics of an equally 
gated DD in the subgap region in Fig. 7 (b-c). In (b) the 
spacing of the 1-particle energy levels \2b\ < 2|A|, hence, 
the excited state can be observed in the current (arrow 
in Fig. 7). In (c) we increase the tunneling strength 
between the two dots 2\b\ > 2|A| and the excited state 
line is vanishing, as explained in Fig. 10. As in the case 
for 2\b\ < 2|A| the excited level is in resonance with the 
St+ transition, however, due to the larger level spacing, 
the ground state transition opens and current is blocked. 



2. Independently gated dots 

In the last paragraph we considered a DD with both 
dots coupled to the same gate electrode. In most ex- 
periments, however, it is more convenient to couple the 
dots independently, which leads to a 'honeycomb' shaped 
current voltage characteristics 37 . For symmetric on-site 
energies and Coulomb repulsion it is possible to diagonal- 
ize the DD Hamiltonian of Eq. (3) analytically. Gating 
the dots independently destroys this symmetry, an ana- 
lytically diagonalization is not possible, and one has to 
use numerical methods. We plot the current as a func- 
tion of the detuning A^ = Vg — Vg , and the average of 
the two gate voltages Y> g = (Vg 1 + V g 2 )/2. 

The current voltage characteristic for serial and par- 
allel configuration is depicted for the normal conducting 
case in Fig. 11 (a)-(b) and for the superconducting case 
in Fig. 11 (c)-(d). Comparing both configurations, we 
observe for the serial one a decrease in the current for 
high detuning A p , while in the parallel configuration cur- 
rent can be observed over the entire voltage range. This 



12 



(a) 



-4 -2 2 4 

eZ g [meV] 



60 < 
40 | 

o 

20 




-2 2 
e£ g [meV] 




-2 2 

eE g [meV] 



Figure 11. (Color online) (a)-(b) Current voltage character- 
istics of a DD coupled to normal conducting leads in serial 
(a) and in parallel (b) configuration. We fixed the bias volt- 
age to eVb = 0.3 meV. (c)-(d) Current voltage characteristics 
of a DD coupled to superconducting leads in serial (c) and 
in parallel (d) configuration. We fixed the bias voltage to 
eVb = 0.3 me V + 2|A| in order to obtain the same conditions 
as for the normal conducting case in (a)-(b). 



difference is a consequence of the geometry of the set- 
up as the DD system remains unchanged. An increase of 
the detuning leads to a localization of the systems ground 
state at site 1 and transitions through site 2 are blocked. 
Since in serial configuration the right lead is only coupled 
to site 2, the localization of the wave function at site 1 
leads to a decrease in the current. In parallel configura- 
tion, however, both sites are coupled to both leads and 




eVb + 2|A| 
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Figure 12. (Color online) (a) E g -N diagram of the 0-1- 
particle transition for eVb/2 > |A|. In the 1-particle spec- 
trum we plotted two situations which mark the borders of 
the current step. The dashed levels mark the left border (for 
small S^) where the 1-particle levels lie above the 0-particle 
energy level. If the distance AE g < eVb/2 — |A| current can 
flow through S+ and D- transitions. By lowering eE g the 1- 
particle energy levels move down in the E g -N diagram, while 
the transitions remain open. The solid lines mark the right 
border of the current steps, as for levels lying below the solid 
line the D- transition is closed and current is blocked. Thus, 
the width of the current steps in the current voltage char- 
acteristics is: eAY> g = eVb — 2|A|. (b) E g -N diagram of 
the 0- 1-particle transition involving thermal transitions. For 
the same arguments as in (a), the distance between two ther- 
mal lines in the current voltage characteristics is equal to 
eAE«7 = eVb + 2|A|. 




1 N 

Figure 13. (Color online) E g -N diagram for the 0- 1-particle 
transition. Transitions between the two 1-particle levels 
(dashed lines) and the 0-particle ground state are allowed 
through the thermal St+ transition. Increasing the gate volt- 
age the levels move down in energy (solid lines) and the ex- 
cited state transition can be observed when the excited level 
aligns with the St+ transition. Hence, the distance of two 
neighboring thermal transitions is equal to the level spacing. 



the ground state transition is always open. 

The left and right border of the current steps are given 
by the source and drain lines, respectively. They follow, 
in complete analogy to the simplest case, from energy 
conservation. In Fig. 12 (a) we show the E g -N diagram 
for the to 1-particle transition illustrating two limits: 
the ground states are (i) in resonance with the S+ tran- 
sition (dashed line) and (ii) in resonance with the D- 
transition (solid line), describing the left and right bor- 
ders of the current step in Fig. 11 (c-d). Starting at the 
S+ resonance, the energy levels of the 1-particle spectrum 
are moving down in energy by increasing the average gate 
voltage Y>g. Both transitions (S+ and D-) remain open as 
long as the ground state lies in the blue (shaded) region. 
If the ground state lies below the solid line, the D- transi- 
tion is closed and current is blocked. Hence, the width of 



the current steps in the current voltage characteristics in 
Fig. 11 (c-d) is equal to the size of the blue (shaded) re- 
gion in Fig. 12 (a), namely eATg = eV& — 2|A|. The same 
arguments hold for the distance of two corresponding 
thermal transitions, as illustrated in Fig. 12 (b) the dis- 
tance of two thermal lines is equal to eAT g = eV& + 2| A|. 

As we can see in Fig. 11 there exists a one to one 
correspondence of the transport conditions of the nor- 
mal conducting to the superconducting case which leads 
to the same shape of the current voltage characteristics 
if ksT <C |A|. Increasing the bias voltage by 2|A| com- 
pared to the normal conducting case eV£ c = eF 6 NC + 2|A| 
leads to the same transport conditions. Although the 
shape of the current steps in Figs. 11 (a-b) and 11 (c- 
d) look the same, they differ at the edges of the current 
steps, as in the superconducting case the sharp peaks 
of the quasiparticle density of states are reflected in the 
current. 



3. Thermal effects 

We have seen that the shape of the stability diagram 
can be explained using energy conservation, in complete 
analogy to the simplest case. In this subsection we like 
to discuss the case for small bias voltages < |A|, 

where current can flow due to thermally excited quasipar- 
ticles exclusively. As already observed above, thermally 
excited quasiparticles do not produce steps in the current 
voltage characteristics rather they appear as small peaks. 
This can be used to resolve transitions through excited 
system states whose energy difference to the ground state 
is less than 2|A|. By detuning the gate voltages of the 
two sites of the DD we can change the level spacing of the 
systems eigenenergies; hence, the excited states are only 
observed in a certain detuning range. To analyze tran- 
sitions through excited system states, c.f. Fig. 14, we 
choose the parallel configuration to rule out the geomet- 
rical effect also leading to a decrease of the current for 
high detuning. If a line corresponding to an excited state 
disappears for higher detuning A gi we conclude that the 
energy difference to its ground state is larger than 2|A|. 
In Fig. 15 we plotted the energy differences of the excited 
states with respect to their ground state for different val- 
ues of the detuning A gi which are marked as red lines in 
Fig. 14. Counting the number of levels lying under the 
red line in Fig. 15 gives information about the number 
of visible excited lines. For instance, consider the case 
of A g = in Fig. 15. Following the red line from small 
to high Tig in Fig. 14, we cross the 0-1 particle transi- 
tions and observe three lines: two corresponding to the 
ground state, and one line in between corresponds to a 
transition through the 1-particle excited state. The dis- 
tance between the leftmost ground state transition line 
and the excited line determines the level spacing of the 
one particle spectrum, see 13. In the 2-particle spectrum 
the energy difference of one excited state lies under the 
red line. Hence we should see two lines coming from ex- 
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Figure 14. (Color online) Current voltage characteristics of 
a DD in parallel configuration for bias V& < 2|A|/e. Since 
the bias voltage is not high enough current can flow only due 
to thermally excited quasiparticles. The red lines correspond 
to Fig. 15 where the energy differences of the excited states 
with respect to their ground state are plotted as a function 
of particle number. The number of visible excited states is 
proportional to the number of energy differences which are 
smaller than 2|A| (red line in Fig. 15). 



cited system states, namely the transition between the 
1-particle ground state and the 2-particle excited state, 
and transitions between the 2-particle ground state and 
the 1-particle excited state. Along the horizontal cut at 
A^ = 2 in Fig. 14, excited states can only be observed 
for the 1-2 particle and the 2-3 particle transition. This 
is in agreement with Fig. 15, where only in the 2 parti- 
cle subspace energy differences lie under the threshold of 
2|A|. For higher detuning, e.g. A^ = 4, no excited states 
can be seen, as the detuning increases the level spacing, 
and all energy differences are larger than 2|A| Fig. 15. 



C. The N-QD-S junction 

We like to close this article by investigating a so called 
N-QD-S hybrid system, where a quantum dot system is 
coupled to a normal and to a superconducting lead, giv- 
ing a possible explanation for the subgap features in Ref. 
1. In the experiment of Ref. 1 a carbon nanotube was 
contacted to two normal conducting leads and to a su- 
perconducting finger in between. The differential conduc- 
tance between the superconducting finger and a normal 
lead is measured, realizing a N-QD-S hybrid system. It 
is possible to apply a bias voltage across the entire tube 
as well as between the superconductor and a normal con- 
ducting lead. The stability diagram in Fig. 2 (a) in Ref. 
1, with no bias applied over the entire tube, reveals the 
typical Coulomb diamond pattern resulting from quasi- 
particle tunneling with no subgap features. By applying 
a bias voltage Vsd over the entire tube, the gap in the sta- 
bility diagram gets smaller with respect to the unbiased 
case and conductance lines can be seen in the Coulomb 
blockade region, c.f. Fig. 3 (a) of Ref. 1. The reduction 
of the gap in the stability diagram is proportional to the 
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Figure 15. (Color online) Plot of the energy differences of the 
excited system states with respect to their ground state as a 
function of particle number. If the energy difference is smaller 
than 2 1 A|, transitions through these excited states can be seen 
in the current voltage characteristics. The threshold of 2|A| 
is marked as a red horizontal line. We depicted the plots for 
three situations differing in the detuning A g . The three cases 
are marked as horizontal lines in Fig. 14. 



applied bias voltage of approximately bVsd ~ I A 1/2, and 
is related to an effective reduction of the superconducting 
gap. For a smaller gap quasiparticles can get thermally 
excited across the gap leading to subgap transport in 
complete analogy to the S-QD-S case discussed above. 

We can model the N-QD-S system by setting |As| = 
for the normal conducting lead (source) in the mas- 
ter equation; the drain contact remains superconduct- 
ing | Ad I = |A|. Hence, the transport conditions change 
slightly and can be summarized in the scheme of Fig. 
18. In Fig. 16 we schematically sketched the expected 
transition lines for a N-QD-S hybrid structure. In Fig. 
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Figure 16. (Color online) Sketch of the transition line of a 
QD coupled to a normal conducting (source) and a supercon- 
ducting lead (drain). The difference to the S-QD-S system 
is that only the drain lines split due to the superconducting 
gap, the S+ and S- lines are described by the same equation. 
In this case a gap equal to |A| is opening, and the triangles 
are shifted apart. Thermal lines can be observed only for the 
drain. 



19 we analyzed the two most important cases, marked as 
points (a) and (b) in Fig. 16. Point (a) shows a paradox- 
ical situation as the particle number of the system seems 
to be increased only at the drain contact, which would 
lead to a negative current at positive bias. However, if 
the two contacts have the same temperature, the thermal 
broadening of the S+ line gives a small contribution in 
the transition rates (dashed green arrow in Fig. 19 (a)) 
making the current positive. The situation in (b) shows 
again the system being in thermal equilibrium with the 
source contact. 

We can see that the lines with negative slope (drain 
lines) give a finite current in the Coulomb blockade re- 
gion as observed in Fig. 3 (b) in the experiments. Thus, 
we claim that the subgap features observed in the exper- 
iments possibly are transitions involving thermally ex- 
cited quasiparticles which are allowed due to the reduc- 
tion of the superconducting gap. This argument is sup- 
ported by the observation that for diamonds where the 
gap has the same size as before (edges of the stability 
diagram), no subgap lines can be observed. In Fig. 17 
we show two dl/dV— characteristic of a N-QD-S system 
corresponding to different superconducting gaps with the 
same temperature (fc#T = O.lmeV) in both cases. In (b) 
the superconducting gap (|A| = 0.3meV) is only half of 
the gap in (a) (fc#T = 0.6meV). By reducing the gap, the 
temperature becomes large enough to excite quasiparti- 
cles across the gap, leading to conductance peaks in the 
Coulomb blockade region, as observed in the experiments 
However, a more complex modeling of the multi-terminal 
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Figure 17. (Color online) Differential Conductance of a SD 
coupled to a normal conducting (source) and to a super- 
conducting lead (drain) (N-QD-S system). The coupling 
to the lead is eT — 0.01 meV. (a) Superconducting gap of 
|A| = 0.6 meV and temperature ksT = O.lmeV. No thermal 
lines in the subgap region are visible, (b) The same tem- 
perature ksT = O.lmeV, but for smaller gap |A| = 0.3meV; 
quasiparticles get thermally excited across the gap leading to 
transport in the Coulomb blockade region. 



eV b < |A| 




Figure 18. (Color online) Visualization of the transport condi- 
tions for a N-QD-S system with eVb/2 < |A|, where the source 
is a normal and the drain a superconducting lead. They fol- 
low from Eqs. (58)- (65) by setting |A| = in the equations 
corresponding to the source lead. 



system is required to understand the experiments in all 
details. 



VI. CONCLUSION 

In this work we developed a transport theory for nanos- 
tructures coupled to superconducting leads up to second 
order in the tunneling Hamiltonian. We used the Bo- 




(b) 



mqL-eV b /2 + |A| 



1 N 



s+ 



eV b /2 



1 N 



Figure 19. E g -N diagrams corresponding to points (a) and 
(b) of Fig. 16. (a) We see a positive current in the subgap 
region, which comes only due to the thermal smearing of the 
S+ transition, (b) The line connecting the S+ and the S- 
transition line in the Coulomb blockade region the system is 
in thermal equilibrium with the source contact. 



goliubov transformation to describe the electrons in the 
superconductors as Cooper pairs and Bogoliubov quasi- 
particle excitations, whereby we modified the Bogoliubov 
transformation in a number conserving way 25 ' 26 , intro- 
ducing Cooper pair creation and annihilation operators 
explicitly. We showed the predictions of the theory on 
two examples, the well known single level quantum dot, 
and the double quantum dot. The characteristic gap in 
the Coulomb diamonds, proportional to the supercon- 
ducting gap, as well as negative differential conductance 
was observed in both cases. Further, we considered the 
double quantum dot in serial as well as in parallel con- 
figuration, see Fig. 1, coupling the dots to the same as 
well as to two separate gate electrodes. 

We systematically analyzed the stability diagrams, ex- 
tending the scheme of Ref. 36 for superconducting leads. 
We found that transport through excited system states 
occurs even for low bias voltages using thermally excited 
quasiparticles, leading to zero bias peaks in the conduc- 
tance. Transitions through excited states can be observed 
if transitions through the ground state are energetically 
not allowed, namely if the distance between the energy 
levels of the excited state and the ground state is smaller 
than 2|A|. This effect can be seen in the the current 
voltage characteristics of an independently gated double 
quantum dot in parallel configuration without tuning pa- 
rameters of the system, since the level spacing changes 
with the detuning A^ of the gate voltages. Hence the 
excited states can be seen only in certain detuning win- 
dows. Finally, we analysed the case where a quantum 
dot is coupled to a normal and a superconducting lead, 
giving a possible explanation for the subgap features of 
Ref. 1 in terms of transport involving thermally excited 
quasiparticles. 

We conclude with the observation that thermally ex- 
cited quasiparticles can lead to a finite current in the 
Coulomb blockade region. Besides the well known ther- 
mal transitions through the ground states, transitions 
through excited system states must be taken into account 
as they are an additional source of zero bias peaks in the 
conductance. For a better comparison with experiments 
the theory can be used to investigate more realistic sys- 
tems like carbon nanotube quantum double dots. Specif- 
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ically, the current voltage spectroscopy in the low bias 
regime can be used to learn something about the spec- 
trum of the set-up. In order to observe the Josephson 
effect and Andreev reflections, the theory must be ex- 
tended to higher order perturbation theory 20 ' 22 . 
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Appendix A: The BCS ground state and the Cooper 
pair operators 

1. Cooper pair states 

In Sect. Ill we claimed that a state \n) contains 2n 
electrons: 



of the Bogoliubov transformation. Combining Eq. (A4) 
with Eq. (20) we find that 



[7V,S f ] =2S f , 
[N,S] = -2S. 



(A7) 



In Eq. (21) we have already seen that S \n) = \n + 1); 
so we still have to proof the action of S |n), using Eqs. 
(A7) and (Al): 

NS\n) = (-2S + SA0 \n) = 2(n - 1) S \n) =N\n-l). 

(A8) 

Hence we proofed that S \n) = \n — l). Acting on a test 
function we finally find: 

S S f |GS) = E b n S \n + 1) = 1 |GS) = S f S |GS) , (A9) 

n 

and consequently 



[S ,S] |GS> =0. 



(A10) 



N \n) = 2n \n) . 



(Al) 



In order to proof this statement we need the following 



commutators for = cj^c^ fc ^: 



[Ri,Rg] = *fcg(4t 6 fct + 6 - 



H c_ H -l), 



^[2B} k + B} k N)\n) 



(A2) 



(A3) 



[N,R\}=2-k\. (A4) 
In the following we proof Eq. (Al) by induction: 
Basis: n = 1 

JV|l> = Jvl£^Rl|0) = 2|l). (A5) 



Inductive step: 

iV>+l>=A/l£^Ri,n> 

= 2(n+l)|n+l>. □ 



(A6) 



2. Properties of Cooper pair operators 

In this paragraph we proof the properties of the Cooper 
pair operators which are necessary to fulfill the unitarity 



Appendix B: Thermodynamic properties of the 
superconducting leads 

In Sect. Ill A we introduced with Eq. (27) the Fermi 
function: 



f + (E) = Tv B (jl% <T p B ). 



(Bl) 



The situation is similar to that of a non-interacting Fermi 
gas, but with the BCS-ground state playing the role of 
the physical vacuum state. Therefore, the derivation is 
standard, see e.g. 38 . Moreover, we have to calculate ther- 
mal expectation values containing an additional Cooper 
pair operator: 

Tr B ( PBTLTfeaSj =^({n}\ fella IW) 

V J in} 

= E \W})({n f }\ S f |{n}> 

{n},{n'} 

= E (WI^7L7^IK})^(GS|S f |GS) 

{n},{n'} 

= Tr B (p B jl% a ) (GS|S f |GS). 



(B2) 

where the trace is over the many-body states defined in 
Eq. (28). The remaining task is to calculate the BCS- 
ground state expectation value of the Cooper pair oper- 
ator: 

(GS|S f |GS) = E 

b n b m (n\ S f \m) 

(B3) 

= ^2 ^n+l^n- 
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Inserting the definition of the b n , we obtain 



(GS| |GS) = Af 2 ^2 



(01) 



2n+l 



n!(n+ 1)! 



Ji(2oi) 
/ (2oi)' 



(B4) 



where we used that Af 2 = ^2 n = io(2ai); io an d 
Ji are modified Bessel functions of the first kind 39 . To 
calculate the actual value for the expectation value, we 
have to evaluate a±: 



k 

pN 



\Uk\ 



/CO 
dEG(\E\ - \A 
-CO 

f 

2p N \A 



\E\ 



\Uk\ 



^E* - [Ap M 2 
/|A| V^ 2 -|A| 2 |£7| + VE 2 ~ |A|2 



, / 2a; 3 
ax ; 

i \V^i 



2x z 



= 2p N \A 
_ 4 Vmk F 

~~ 3 7T 2 /L 2 



X 2 - 1 - -or 3 



IAI. 



(B5) 



The elementary integrals can be found e.g. in Ref. 40. 
In the last step of Eq. (B5) we calculated the limit x — >> 
oo by Taylor expanding the roots. We see that a\ is 
proportional to W, hence we find in the thermodynamic 
limit: 



(GS| S f |GS) = ^4 ^ (B6) 

io(2ai) 



using the asymptotic expansion of the modified Bessel 
functions 41 : 



In(x) 



Finally we find that 



x i 

7=[l + 0(-)]. 



(B7) 



(B8) 



Appendix C: Rates 
1. Normal rates 



In the stationary limit, r — » oo, the normal rates read: 

2 

N—^N-\-l " ' " 



(ra|d ao . |m) (m'ld^ |n) / 
Jo 



dt 2 e^ En '™' t2 



(CI) 



(n| dL |m) (m'| d aV |n') / dt 2 e* B »'™'* 2 
Jo 

K fc r/-(^*)e-* (B " 1+ '"' ) * a + |^ fe | 2 / + (i3 t , fe )e + ^^-^)* 2 



(C2) 



2. Anomalous rates 



The anomalous rates are 



(n| dL |m) (m'| d^ |n'> f dt 2 e i B ^ e "i 2 ^ 

/ + (£:, 7fc )e + s ( ' E '" t+A1 '' ) * 2 - f-(E nk )e-i {E ^- fl '' )t2 , 

(C3) 

(^nmm'n^^^" 1 = T !™Q) S Sgn (a) 



f + (E vk )e + ^ (E ^ k ~^ )t2 - f-(E nk )e~^ Er ' k+il '> )t2 

(C4) 

From Eqs. (C3) and (C4) we can see that they contain 
an integral of the form: 



dEF(E) 



f^\E)e^ qE+uj)t2 



(C5) 



0, 



where F(—E) = F(£") and p,q = ±1. The anomalous 
rates with the superscript '— ' contain the same structure, 
hence all anomalous rates are equal to zero. 
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3. Renormalization of the rates 

In the lowest order approximation we find rates which 
are proportional to the BCS-density of states leading to 
divergences at the gap edges. We can can renormalize 
the rates by introducing a finite lifetime (j/h) -1 in the 
exponents of Eqs. (CI) and Eqs. (C2). Since we are 
neglecting coherences the imaginary parts of the rates 
do not contribute to the dynamics of the system. For 
example consider the integral appearing in Eq. (CI): 

dE J™ dt 2 e^ E+UJ+i ^f + (E)D(E)j 



dE 



(C6) 



f+(E)D(E), 



where we introduced uo = E n > m > + /i^. Generalizing the 
integral for the cases (TV — >• N ± 1) it reads 



/CO rOO 
dE L(E,lj) f ± (E) D(E) = ft / dEF(E), 
-co J —oo 



where 



L(E,u) 



7 



(£ + c^) 2 + 7 2 



(C7) 



(C8) 



describes the Lorentzian and F(E) = 
L(E,ou) f ± (E) D(E). We can solve the integral of 
Eq. (C7) using residue calculus hence. To this extend 
we analyze the singularities of the integrand and the 
area in which the integrand is analytic. The Lorentzian 
L(E,cj) has poles at 



E = —uj =p Z7, 



with the corresponding residues: 



ResE= 



-u;=F7 



L(E) 



±i 



(C9) 



(CIO) 



The poles of the Fermi function f ± {E) are purely imag- 
inary and equally distributed along the imaginary axis: 



E= — (2n + l) neZ, 



with the residues 
Res 



B=3f (2n+l) / ± (-^') 



(Cll) 



(C12) 



The square roots in the BCS-density of states D(E) have 
branch cuts along the real axis. In Fig. 20 we sketched 
the contour in the complex plane which is slightly shifted 
away from the real axis with e = 1/R. In the limit R — >• 
oo the integral along the semicircle vanishes and we are 
left with: 



lim 



dxF(x + ie) = 2m V] Res z=a F(z). (C13) 

-R 

lm[zY 




R Re(z) 



Figure 20. Contour in the complex plane used to integrate 
Eq. (C7). 



In the limit R —> oo Eq. (C13) is mapped back into the 
real integral of Eq. (C7), and we find: 

/CO 
dE L(E) f ± (E)D(E) 

~°7 (C14) 
=irhRe ( f + (-uj + ij) D(E - uo + ij) J . 
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